#include "Metaball.h"
#include <algorithm>
#include "Polynom.h"



Metaball::Metaball(std::vector<VECTOR> acenters, std::vector<SCALAR> aradius, double am_thershold)
{
	centers = acenters;
	radius = aradius;
	fieldsPolynoms.resize(centers.size());
	usedPolys.resize(centers.size());
	m_thershold = am_thershold;
}

//take the polynom of the function of the field density for one sphere in the metaball and one ray
void Metaball::getPolysFieldDenisty(const CLine & line, int metaballID)
{
	VECTOR center = centers[metaballID];
	SCALAR invRadius = 1 / radius[metaballID];
	VECTOR rayDir = line.dir;
	VECTOR D = line.loc - center;
	SCALAR invSquareradius = invRadius*invRadius;
	SCALAR invQuadradius = invSquareradius*invSquareradius;
	double coef[3];
	//determine dist function depending of t
	coef[0] = D.x * D.x + D.y * D.y + D.z * D.z ;
	coef[1] =    2 * (D.x * rayDir.x + D.y * rayDir.y + D.z * rayDir.z);
	coef[2] = 1;
	//determine field function depending of t
	//***********************************//
	//I(t) = dist^2/r^4 - 2*dist/r^2 +1**//
	//***********************************//
	Polinomio p2 = Polinomio(coef , 2);
	Polinomio p = Polinomio(coef , 2);
	p2 * p2;
	p2*invQuadradius;
	double b = invSquareradius * 2;
	p*b;
	Polinomio field = p2-p;
	field.sumCoef(0,1.0f);
	
	//asingn the polynom to metaball
	fieldsPolynoms[metaballID] = field;
}

Intersections Metaball::getMetaballsIntersections(const CLine & line)
{
	
	Intersections intersections = Intersections();
	bool intersects;
	double t0;
	double t1;
	for(unsigned i = 0; i < centers.size(); i ++){
		CSphere  currentMetaball = CSphere(radius[i]);
		currentMetaball.setLocation(centers[i]);
		intersects = currentMetaball.hits(line, t0, t1);
		if(intersects){
			intersections.push_back(Intersection(t0, true, i));
			intersections.push_back(Intersection(t1, false, i));
		}
	}	
	return intersections;
}

bool sortIntersections(const Intersection& i1, const Intersection& i2)
{
	if(i1.t==i2.t)
		return i1.in;
	else
		return i1.t < i2.t;
}


bool Metaball::hits (const CLine &line, SCALAR &t_hit)
{
	//get the ray intersections with the spheres 
	Intersections intersections = getMetaballsIntersections(line);
	
	//there isn't intersection
	if(intersections.size()<2)
		return false;
	
	//sort intersections
	std::sort(intersections.begin(), intersections.end(), sortIntersections);

	//init the polynom and other parameters;
	SCALAR tMin;
	SCALAR tMax;
	double initialCoefs[5] = {-m_thershold,0,0,0,0};
	Polinomio fieldPoly = Polinomio(initialCoefs,4);
		
	//change the fieldPolynom with the intersections and search the roots
	for(unsigned intIt = 0; intIt < intersections.size()-1; intIt++)
	{
		Intersection currentIntersection = intersections[intIt];
		Intersection nextIntersection = intersections[intIt+1];
		tMin = currentIntersection.t;
		tMax = nextIntersection.t;
		int metaballID = currentIntersection.metaballId;
	
		if(currentIntersection.in==true)
		{
			getPolysFieldDenisty(line, metaballID);
			fieldPoly = fieldPoly + fieldsPolynoms[metaballID];
		}
		else
			{
			fieldPoly = fieldPoly - fieldsPolynoms[metaballID];
			}
		std::vector<double> coefs;
		coefs.resize(5);
		coefs[0]=fieldPoly.Coef(0);
		coefs[1]=fieldPoly.Coef(1);
		coefs[2]=fieldPoly.Coef(2);
		coefs[3]=fieldPoly.Coef(3);
		coefs[4]=fieldPoly.Coef(4);
		//search the roots of the polynom
		std::vector<double> roots = test_polynomial(coefs);
		std::sort(roots.begin(),roots.end());
		for(unsigned rootIt = 0; rootIt < roots.size(); rootIt++)
		{
			double root = roots[rootIt];
			if(root > tMin+SMALL_AMOUNT && root < tMax-SMALL_AMOUNT){//is the root in the interval tmin, tmax?
				t_hit = root;//intersect!
				tHit = root;
				return true;
			}
		}
	}
	return false;
}

VECTOR Metaball::getNormal(const VECTOR &hit_loc) {
	VECTOR gradient = VECTOR(0 , 0 , 0);
	//calculate the gradient of the field function for the hit_loc point
	for(unsigned i = 0; i < centers.size(); i++){
			VECTOR current_normal = hit_loc - centers[i];
			SCALAR InvSquareNormalDotNormal = current_normal.dot(current_normal);
			InvSquareNormalDotNormal = 1 / (InvSquareNormalDotNormal * InvSquareNormalDotNormal);
			current_normal = 2 * InvSquareNormalDotNormal * (centers[i] - hit_loc);
			gradient += current_normal;
		
	}
	gradient = gradient.normalize();

	if(material->bumpMap)
	{
		 material->bumpNormal(hit_loc,gradient);
	}
	return gradient;
}

